      SUBROUTINE WNNLS(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE,    WNN   10
     * IWORK, WORK)
C
C     DIMENSION W(MDW,N+1),PRGOPT(*),X(N),IWORK(M+N),WORK(M+5*N)
C
C     ABSTRACT
C
C     THIS SUBPROGRAM SOLVES A LINEARLY CONSTRAINED LEAST SQUARES
C     PROBLEM.  SUPPOSE THERE ARE GIVEN MATRICES E AND A OF
C     RESPECTIVE DIMENSIONS ME BY N AND MA BY N, AND VECTORS F
C     AND B OF RESPECTIVE LENGTHS ME AND MA.  THIS SUBROUTINE
C     SOLVES THE PROBLEM
C
C               EX = F, (EQUATIONS TO BE EXACTLY SATISFIED)
C
C               AX = B, (EQUATIONS TO BE APPROXIMATELY SATISFIED,
C                        IN THE LEAST SQUARES SENSE)
C
C               SUBJECT TO COMPONENTS L+1,...,N NONNEGATIVE
C
C     ANY VALUES ME.GE.0, MA.GE.0 AND 0.LE. L .LE.N ARE PERMITTED.
C
C     THE PROBLEM IS REPOSED AS PROBLEM WNNLS
C
C               (WT*E)X = (WT*F)
C               (   A)    (   B), (LEAST SQUARES)
C               SUBJECT TO COMPONENTS L+1,...,N NONNEGATIVE.
C
C     THE SUBPROGRAM CHOOSES THE HEAVY WEIGHT (OR PENALTY PARAMETER) WT.
C
C     THE PARAMETERS FOR WNNLS ARE
C
C     INPUT..
C
C     W(*,*),MDW,  THE ARRAY W(*,*) IS DOUBLE SUBSCRIPTED WITH FIRST
C     ME,MA,N,L    DIMENSIONING PARAMETER EQUAL TO MDW.  FOR THIS
C                  DISCUSSION LET US CALL M = ME + MA.  THEN MDW
C                  MUST SATISFY MDW.GE.M.  THE CONDITION MDW.LT.M
C                  IS AN ERROR.
C
C                  THE ARRAY W(*,*) CONTAINS THE MATRICES AND VECTORS
C
C                       (E  F)
C                       (A  B)
C
C                  IN ROWS AND COLUMNS 1,...,M AND 1,...,N+1
C                  RESPECTIVELY.  COLUMNS 1,...,L CORRESPOND TO
C                  UNCONSTRAINED VARIABLES X(1),...,X(L).  THE
C                  REMAINING VARIABLES ARE CONSTRAINED TO BE
C                  NONNEGATIVE.  THE CONDITION L.LT.0 .OR. L.GT.N IS
C                  AN ERROR.
C
C     PRGOPT(*)    THIS ARRAY IS THE OPTION VECTOR.
C                  IF THE USER IS SATISFIED WITH THE NOMINAL
C                  SUBPROGRAM FEATURES SET
C
C                  PRGOPT(1)=1 (OR PRGOPT(1)=1.0)
C
C                  OTHERWISE PRGOPT(*) IS A LINKED LIST CONSISTING OF
C                  GROUPS OF DATA OF THE FOLLOWING FORM
C
C                  LINK
C                  KEY
C                  DATA SET
C
C                  THE PARAMETERS LINK AND KEY ARE EACH ONE WORD.
C                  THE DATA SET CAN BE COMPRISED OF SEVERAL WORDS.
C                  THE NUMBER OF ITEMS DEPENDS ON THE VALUE OF KEY.
C                  THE VALUE OF LINK POINTS TO THE FIRST
C                  ENTRY OF THE NEXT GROUP OF DATA WITHIN
C                  PRGOPT(*).  THE EXCEPTION IS WHEN THERE ARE
C                  NO MORE OPTIONS TO CHANGE.  IN THAT
C                  CASE LINK=1 AND THE VALUES KEY AND DATA SET
C                  ARE NOT REFERENCED. THE GENERAL LAYOUT OF
C                  PRGOPT(*) IS AS FOLLOWS.
C
C               ...PRGOPT(1)=LINK1 (LINK TO FIRST ENTRY OF NEXT GROUP)
C               .  PRGOPT(2)=KEY1 (KEY TO THE OPTION CHANGE)
C               .  PRGOPT(3)=DATA VALUE (DATA VALUE FOR THIS CHANGE)
C               .       .
C               .       .
C               .       .
C               ...PRGOPT(LINK1)=LINK2 (LINK TO THE FIRST ENTRY OF
C               .                       NEXT GROUP)
C               .  PRGOPT(LINK1+1)=KEY2 (KEY TO THE OPTION CHANGE)
C               .  PRGOPT(LINK1+2)=DATA VALUE
C               ...     .
C               .       .
C               .       .
C               ...PRGOPT(LINK)=1 (NO MORE OPTIONS TO CHANGE)
C
C                  VALUES OF LINK THAT ARE NONPOSITIVE ARE ERRORS.
C                  A VALUE OF LINK.GT.NLINK=100000 IS ALSO AN ERROR.
C                  THIS HELPS PREVENT USING INVALID BUT POSITIVE
C                  VALUES OF LINK THAT WILL PROBABLY EXTEND
C                  BEYOND THE PROGRAM LIMITS OF PRGOPT(*).
C                  UNRECOGNIZED VALUES OF KEY ARE IGNORED.  THE
C                  ORDER OF THE OPTIONS IS ARBITRARY AND ANY NUMBER
C                  OF OPTIONS CAN BE CHANGED WITH THE FOLLOWING
C                  RESTRICTION.  TO PREVENT CYCLING IN THE
C                  PROCESSING OF THE OPTION ARRAY A COUNT OF THE
C                  NUMBER OF OPTIONS CHANGED IS MAINTAINED.
C                  WHENEVER THIS COUNT EXCEEDS NOPT=1000 AN ERROR
C                  MESSAGE IS PRINTED AND THE SUBPROGRAM RETURNS.
C
C                  OPTIONS..
C
C                  KEY=6
C                         SCALE THE NONZERO COLUMNS OF THE
C                  ENTIRE DATA MATRIX
C                  (E)
C                  (A)
C                  TO HAVE LENGTH ONE.  THE DATA SET FOR
C                  THIS OPTION IS A SINGLE VALUE.  IT MUST
C                  BE NONZERO IF UNIT LENGTH COLUMN SCALING IS
C                  DESIRED.
C
C                  KEY=7
C                         SCALE COLUMNS OF THE ENTIRE DATA MATRIX
C                  (E)
C                  (A)
C                  WITH A USER-PROVIDED DIAGONAL MATRIX.
C                  THE DATA SET FOR THIS OPTION CONSISTS
C                  OF THE N DIAGONAL SCALING FACTORS, ONE FOR
C                  EACH MATRIX COLUMN.
C
C                  KEY=8
C                         CHANGE THE RANK DETERMINATION TOLERANCE FROM
C                  THE NOMINAL VALUE OF DSQRT(EPS).  THIS QUANTITY CAN
C                  BE NO SMALLER THAN EPS, THE ARITHMETIC-
C                  STORAGE PRECISION.  THE QUANTITY USED
C                  HERE IS INTERNALLY RESTRICTED TO BE AT
C                  LEAST EPS.  THE DATA SET FOR THIS OPTION
C                  IS THE NEW TOLERANCE.
C
C                  KEY=9
C                         CHANGE THE BLOW-UP PARAMETER FROM THE
C                  NOMINAL VALUE OF DSQRT(EPS).  THE RECIPROCAL OF
C                  THIS PARAMETER IS USED IN REJECTING SOLUTION
C                  COMPONENTS AS TOO LARGE WHEN A VARIABLE IS
C                  FIRST BROUGHT INTO THE ACTIVE SET.  TOO LARGE
C                  MEANS THAT THE PROPOSED COMPONENT TIMES THE
C                  RECIPROCAL OF THE PARAMETERIS NOT LESS THAN
C                  THE RATIO OF THE NORMS OF THE RIGHT-SIDE
C                  VECTOR AND THE DATA MATRIX.
C                  THIS PARAMETER CAN BE NO SMALLER THAN EPS,
C                  THE ARITHMETIC-STORAGE PRECISION.
C
C                  FOR EXAMPLE, SUPPOSE WE WANT TO PROVIDE
C                  A DIAGONAL MATRIX TO SCALE THE PROBLEM
C                  MATRIX AND CHANGE THE TOLERANCE USED FOR
C                  DETERMINING LINEAR DEPENDENCE OF DROPPED COL
C                  VECTORS.  FOR THESE OPTIONS THE DIMENSIONS OF
C                  PRGOPT(*) MUST BE AT LEAST N+6.  THE FORTRAN
C                  STATEMENTS DEFINING THESE OPTIONS WOULD
C                  BE AS FOLLOWS.
C
C                  PRGOPT(1)=N+3 (LINK TO ENTRY N+3 IN PRGOPT(*))
C                  PRGOPT(2)=7 (USER-PROVIDED SCALING KEY)
C
C                  CALL DCOPY(N,D,1,PRGOPT(3),1) (COPY THE N
C                  SCALING FACTORS FROM A USER ARRAY CALLED D(*)
C                  INTO PRGOPT(3)-PRGOPT(N+2))
C
C                  PRGOPT(N+3)=N+6 (LINK TO ENTRY N+6 OF PRGOPT(*))
C                  PRGOPT(N+4)=8 (LINEAR DEPENDENCE TOLERANCE KEY)
C                  PRGOPT(N+5)=... (NEW VALUE OF THE TOLERANCE)
C
C                  PRGOPT(N+6)=1 (NO MORE OPTIONS TO CHANGE)
C
C     IWORK(1),    THE AMOUNTS OF WORKING STORAGE ACTUALLY ALLOCATED
C     IWORK(2)     FOR THE WORKING ARRAYS WORK(*) AND IWORK(*),
C                  RESPECTIVELY.  THESE QUANTITIES ARE COMPARED WITH
C                  THE ACTUAL AMOUNTS OF STORAGE NEEDED FOR WNNLS( ).
C                  INSUFFICIENT STORAGE ALLOCATED FOR EITHER WORK(*)
C                  OR IWORK(*) IS CONSIDERED AN ERROR.  THIS FEATURE
C                  WAS INCLUDED IN WNNLS( ) BECAUSE MISCALCULATING
C                  THE STORAGE FORMULAS FOR WORK(*) AND IWORK(*)
C                  MIGHT VERY WELL LEAD TO SUBTLE AND HARD-TO-FIND
C                  EXECUTION ERRORS.
C
C                  THE LENGTH OF WORK(*) MUST BE AT LEAST
C
C                  LW = ME+MA+5*N
C                  THIS TEST WILL NOT BE MADE IF IWORK(1).LE.0.
C
C                  THE LENGTH OF IWORK(*) MUST BE AT LEAST
C
C                  LIW = ME+MA+N
C                  THIS TEST WILL NOT BE MADE IF IWORK(2).LE.0.
C
C     OUTPUT..
C
C     X(*)         AN ARRAY DIMENSIONED AT LEAST N, WHICH WILL
C                  CONTAIN THE N COMPONENTS OF THE SOLUTION VECTOR
C                  ON OUTPUT.
C
C     RNORM        THE RESIDUAL NORM OF THE SOLUTION.  THE VALUE OF
C                  RNORM CONTAINS THE RESIDUAL VECTOR LENGTH OF THE
C                  EQUALITY CONSTRAINTS AND LEAST SQUARES EQUATIONS.
C
C     MODE         THE VALUE OF MODE INDICATES THE SUCCESS OR FAILURE
C                  OF THE SUBPROGRAM.
C
C                  MODE = 0  SUBPROGRAM COMPLETED SUCCESSFULLY.
C
C                       = 1  MAX. NUMBER OF ITERATIONS (EQUAL TO
C                            3*(N-L)) EXCEEDED. NEARLY ALL PROBLEMS
C                            SHOULD COMPLETE IN FEWER THAN THIS
C                            NUMBER OF ITERATIONS. AN APPROXIMATE
C                            SOLUTION AND ITS CORRESPONDING RESIDUAL
C                            VECTOR LENGTH ARE IN X(*) AND RNORM.
C
C                       = 2  USAGE ERROR OCCURRED.  THE OFFENDING
C                            CONDITION IS NOTED WITH THE ERROR
C                            PROCESSING SUBPROGRAM, XERROR( ).
C
C     USER-DESIGNATED
C     WORKING ARRAYS..
C
C     WORK(*)      A WORKING ARRAY OF LENGTH AT LEAST
C                  M + 5*N.
C
C     IWORK(*)     AN INTEGER-VALUED WORKING ARRAY OF LENGTH AT LEAST
C                  M+N.
C
C     THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO
C     DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES.
C     USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/.
C     (START AT LINE WITH C++ IN COLS. 1-3.)
C     /REAL (12 BLANKS)/DOUBLE PRECISION/,/, DUMMY/,SNGL(DUMMY)/
C
C     WRITTEN BY KAREN H. HASKELL, SANDIA LABORATORIES,
C     AND R.J. HANSON, SANDIA LABORATORIES.
C     REVISED FEB.25, 1982.
C
C     SUBROUTINES CALLED BY WNNLS( )
C
C++
C     WNLSM         COMPANION SUBROUTINE TO WNNLS( ), WHERE
C                   MOST OF THE COMPUTATION TAKES PLACE.
C
C     XERROR,XERRWV FROM SLATEC ERROR PROCESSING PACKAGE.
C                   THIS IS DOCUMENTED IN SANDIA TECH. REPT.,
C                   SAND78-1189.
C
C     REFERENCES
C
C     1. SOLVING LEAST SQUARES PROBLEMS, BY C.L. LAWSON
C        AND R.J. HANSON.  PRENTICE-HALL, INC. (1974).
C
C     2. BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE, BY
C        C.L. LAWSON, R.J. HANSON, D.R. KINCAID, AND F.T. KROGH.
C        TOMS, V. 5, NO. 3, P. 308.  ALSO AVAILABLE AS
C        SANDIA TECHNICAL REPORT NO. SAND77-0898.
C
C     3. AN ALGORITHM FOR LINEAR LEAST SQUARES WITH EQUALITY
C        AND NONNEGATIVITY CONSTRAINTS, BY K.H. HASKELL AND
C        R.J. HANSON.  AVAILABLE AS SANDIA TECHNICAL REPORT NO.
C        SAND77-0552, AND MATH. PROGRAMMING, VOL. 21, (1981), P. 98-118.
C
C     4. SLATEC COMMON MATH. LIBRARY ERROR HANDLING
C        PACKAGE.  BY R. E. JONES.  AVAILABLE AS SANDIA
C        TECHNICAL REPORT SAND78-1189.
C
      DOUBLE PRECISION  DUMMY, W(MDW,1), PRGOPT(1), X(1), WORK(1), RNORM
      INTEGER IWORK(1)
C
C
      MODE = 0
      IF (MA+ME.LE.0 .OR. N.LE.0) RETURN
      IF (.NOT.(IWORK(1).GT.0)) GO TO 20
      LW = ME + MA + 5*N
      IF (.NOT.(IWORK(1).LT.LW)) GO TO 10
      NERR = 2
      IOPT = 1
      CALL XERRWV(70HWNNLS( ), INSUFFICIENT STORAGE ALLOCATED FOR WORK(*
     *), NEED LW=I1 BELOW, 70, NERR, IOPT, 1, LW, 0, 0, DUMMY, DUMMY)
      MODE = 2
      RETURN
   10 CONTINUE
   20 IF (.NOT.(IWORK(2).GT.0)) GO TO 40
      LIW = ME + MA + N
      IF (.NOT.(IWORK(2).LT.LIW)) GO TO 30
      NERR = 2
      IOPT = 1
      CALL XERRWV(72HWNNLS( ), INSUFFICIENT STORAGE ALLOCATED FOR IWORK(
     **), NEED LIW=I1 BELOW, 72, NERR, IOPT, 1, LIW, 0, 0, DUMMY, DUMMY)
      MODE = 2
      RETURN
   30 CONTINUE
   40 IF (.NOT.(MDW.LT.ME+MA)) GO TO 50
      NERR = 1
      IOPT = 1
      CALL XERROR(44HWNNLS( ), THE VALUE MDW.LT.ME+MA IS AN ERROR, 44,
     * NERR, IOPT)
      MODE = 2
      RETURN
   50 IF (0.LE.L .AND. L.LE.N) GO TO 60
      NERR = 2
      IOPT = 1
      CALL XERROR(39HWNNLS( ), L.LE.0.AND.L.LE.N IS REQUIRED, 39, NERR,
     * IOPT)
      MODE = 2
      RETURN
C
C     THE PURPOSE OF THIS SUBROUTINE IS TO BREAK UP THE ARRAYS
C     WORK(*) AND IWORK(*) INTO SEPARATE WORK ARRAYS
C     REQUIRED BY THE MAIN SUBROUTINE WNLSM( ).
C
   60 L1 = N + 1
      L2 = L1 + N
      L3 = L2 + ME + MA
      L4 = L3 + N
      L5 = L4 + N
C
      CALL WNLSM(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, IWORK,
     * IWORK(L1), WORK(1), WORK(L1), WORK(L2), WORK(L3), WORK(L4),
     * WORK(L5))
      RETURN
      END
